Fit temperature models and predict growing season temperature

Author

Max Lindmark, Jan Ohlberger, Anna Gårdmark

Published

January 15, 2024

Load libraries

pkgs <- c("here","tidyverse", "tidylog", "RColorBrewer", "viridis", "sdmTMB", "sdmTMBextra", "patchwork", "RCurl", "tidylog") 

invisible(lapply(pkgs, library,character.only = T))
here() starts at /Users/maxlindmark/Dropbox/Max work/R/perch-growth
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Attaching package: 'tidylog'


The following objects are masked from 'package:dplyr':

    add_count, add_tally, anti_join, count, distinct, distinct_all,
    distinct_at, distinct_if, filter, filter_all, filter_at, filter_if,
    full_join, group_by, group_by_all, group_by_at, group_by_if,
    inner_join, left_join, mutate, mutate_all, mutate_at, mutate_if,
    relocate, rename, rename_all, rename_at, rename_if, rename_with,
    right_join, sample_frac, sample_n, select, select_all, select_at,
    select_if, semi_join, slice, slice_head, slice_max, slice_min,
    slice_sample, slice_tail, summarise, summarise_all, summarise_at,
    summarise_if, summarize, summarize_all, summarize_at, summarize_if,
    tally, top_frac, top_n, transmute, transmute_all, transmute_at,
    transmute_if, ungroup


The following objects are masked from 'package:tidyr':

    drop_na, fill, gather, pivot_longer, pivot_wider, replace_na,
    spread, uncount


The following object is masked from 'package:stats':

    filter


Loading required package: viridisLite


Attaching package: 'sdmTMBextra'


The following objects are masked from 'package:sdmTMB':

    dharma_residuals, extract_mcmc



Attaching package: 'RCurl'


The following object is masked from 'package:tidyr':

    complete
# devtools::install_github("seananderson/ggsidekick") # not on CRAN 
library(ggsidekick)
theme_set(theme_sleek())

# Set path:
home <- here::here()

Load cache

# qwraps2::lazyload_cache_dir(path = paste0(home, "/R/analyze-data/01-fit-temp-models-predict_cache/html"))

Read data

d <- read_csv(paste0(home, "/data/for-analysis/temp_data_for_fitting.csv"))
Rows: 98616 Columns: 13
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (4): area, source, db, id
dbl  (8): year, temp, yday, month, day, VkDag, stn_nr, section_nr
date (1): date

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
d <- d %>% mutate(area = as.factor(area),
                 source_f = as.factor(source),
                 year_f = as.factor(year)) %>% 
  filter(!area %in% c("VN", "TH")) # VN: no logger data, TH: to short time series
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
filter: removed 5,545 rows (6%), 93,071 rows remaining
# Keep track of the years for which we have cohorts for matching with cohort data
gdat <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/perch-growth/master/data/for-analysis/dat.csv")
Rows: 338460 Columns: 12
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (6): age_ring, area, gear, ID, sex, keep
dbl (6): length_mm, age_bc, age_catch, catch_year, cohort, final_length

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
gdat$area_cohort_age <- as.factor(paste(gdat$area, gdat$cohort, gdat$age_bc))

gdat <- gdat %>%
  group_by(area_cohort_age) %>% 
  filter(n() > 10) %>% 
  filter(age_catch > 3) %>% 
  group_by(area) %>%
  summarise(min = min(cohort),
            max = max(cohort)) %>% 
  arrange(min)
group_by: one grouping variable (area_cohort_age)
filter (grouped): removed 5,190 rows (2%), 333,270 rows remaining
filter (grouped): removed 107,058 rows (32%), 226,212 rows remaining
group_by: one grouping variable (area)
summarise: now 12 rows and 3 columns, ungrouped
d <- left_join(d, gdat, by = "area") %>%
  mutate(area = as.factor(area),
         growth_dat = ifelse(year >= min & year <= max, "Y", "N"))
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (     2)
           > matched rows     93,071
           >                 ========
           > rows total       93,071
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
# Drop data in SI_HA and BT before onset of warming
d <- d %>%
  mutate(discard = "N",
         discard = ifelse(area == "BT" & year <= 1980, "Y", discard),
         discard = ifelse(area == "SI_HA" & year <= 1972, "Y", discard)) %>% 
  filter(discard == "N")
mutate: new variable 'discard' (character) with 2 unique values and 0% NA
filter: removed 1,146 rows (1%), 91,925 rows remaining
# Drop heated areas actually for the full plot
df <- d %>% filter(!area %in% c("BT", "SI_HA"))
filter: removed 24,149 rows (26%), 67,776 rows remaining

Plot data

df %>%
  filter(growth_dat == "Y") %>% 
  distinct(area, year, source) %>%
  ggplot(aes(year, area, color = source)) +
  geom_point(position = position_dodge(width = 0.4), shape = 15) + 
  labs(y = "Site", x = "Year", color = "Source") + 
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "bottom")
filter: removed 24,060 rows (35%), 43,716 rows remaining
distinct: removed 43,092 rows (99%), 624 rows remaining

ggsave(paste0(home, "/figures/supp/temp_sources.pdf"), width = 15, height = 11, units = "cm")

Area-specific models

spec_preds <- list()
res_list <- list()

for(i in unique(d$area)) {
  
  dd <- d %>% filter(area == i)

  if(unique(dd$area) %in% c("BS", "BT", "FB", "FM", "MU", "SI_EK")) { # RA, JM, HO, SI_HA remains
    
    mspec <- sdmTMB(temp ~ 0 + source_f + year_f + s(yday, bs = "cc"), 
                    data = dd,
                    family = student(df = 6),
                    spatial = "off",
                    spatiotemporal = "off",
                    knots = list(yday = c(0.5, 364.5)),
                    control = sdmTMBcontrol(newton_loops = 1)) 
  
  } else {
    
    mspec <- sdmTMB(temp ~ 0 + source_f + year_f + s(yday, bs = "cc"), 
                    data = dd,
                    family = student(df = 10),
                    spatial = "off",
                    spatiotemporal = "off",
                    knots = list(yday = c(0.5, 364.5)),
                    control = sdmTMBcontrol(newton_loops = 1)) 
    
  }
  
  sanity(mspec)

  # Residuals
  mcmc_res_msep <- residuals(mspec, type = "mle-mcmc",
                             mcmc_samples = sdmTMBextra::predict_mle_mcmc(mspec,
                                                                          mcmc_iter = 201,
                                                                          mcmc_warmup = 200))
  
  dd$res <- as.vector(mcmc_res_msep)
    
  # Store residuals
  res_list[[i]] <- dd
  
  print(ggplot(dd, aes(sample = mcmc_res_msep)) +
    stat_qq(size = 0.75, shape = 21, fill = NA) +
    stat_qq_line() +
    labs(y = "Sample Quantiles", x = "Theoretical Quantiles", title = paste("Site = ", i)) + 
    theme(aspect.ratio = 1))
  
  # Predict on new data
  nd_area <- data.frame(expand.grid(yday = seq(min(dd$yday), max(dd$yday), by = 1),
                                    area = unique(dd$area),
                                    source = unique(dd$source_f),
                                    year = unique(dd$year))) %>%
    mutate(source_f = as.factor(source),
           year_f = as.factor(year)) %>% 
    left_join(gdat, by = "area") %>% 
    mutate(area = as.factor(area),
           growth_dat = ifelse(year >= min & year <= max, "Y", "N"))
  
  nd_area$pred <- predict(mspec, newdata = nd_area)$est

  # Save!
  spec_preds[[i]] <- nd_area
  
}
filter: removed 88,464 rows (96%), 3,461 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000395 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 3.95 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 2.30019 seconds (Warm-up)
Chain 1:                0.005236 seconds (Sampling)
Chain 1:                2.30542 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     83,664
           >                 ========
           > rows total       83,664
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 82,953 rows (90%), 8,972 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000996 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 9.96 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 5.6973 seconds (Warm-up)
Chain 1:                0.007769 seconds (Sampling)
Chain 1:                5.70507 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 42 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     46,116
           >                 ========
           > rows total       46,116
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 82,452 rows (90%), 9,473 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001036 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10.36 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 6.26191 seconds (Warm-up)
Chain 1:                0.016907 seconds (Sampling)
Chain 1:                6.27882 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     90,885
           >                 ========
           > rows total       90,885
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 79,957 rows (87%), 11,968 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001176 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 11.76 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 6.92627 seconds (Warm-up)
Chain 1:                0.017853 seconds (Sampling)
Chain 1:                6.94412 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     91,134
           >                 ========
           > rows total       91,134
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 83,763 rows (91%), 8,162 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000863 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 8.63 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 4.34097 seconds (Warm-up)
Chain 1:                0.013235 seconds (Sampling)
Chain 1:                4.3542 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     90,885
           >                 ========
           > rows total       90,885
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 81,791 rows (89%), 10,134 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001076 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10.76 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 6.83227 seconds (Warm-up)
Chain 1:                0.008257 seconds (Sampling)
Chain 1:                6.84053 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     84,660
           >                 ========
           > rows total       84,660
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 81,040 rows (88%), 10,885 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001219 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 12.19 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 7.00604 seconds (Warm-up)
Chain 1:                0.014502 seconds (Sampling)
Chain 1:                7.02054 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     91,134
           >                 ========
           > rows total       91,134
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 87,105 rows (95%), 4,820 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000583 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 5.83 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 3.779 seconds (Warm-up)
Chain 1:                0.007696 seconds (Sampling)
Chain 1:                3.78669 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     91,134
           >                 ========
           > rows total       91,134
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 83,052 rows (90%), 8,873 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001012 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10.12 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 5.75575 seconds (Warm-up)
Chain 1:                0.007634 seconds (Sampling)
Chain 1:                5.76339 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     91,134
           >                 ========
           > rows total       91,134
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 76,748 rows (83%), 15,177 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001591 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 15.91 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 6.97252 seconds (Warm-up)
Chain 1:                0.012057 seconds (Sampling)
Chain 1:                6.98458 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 50 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     54,900
           >                 ========
           > rows total       54,900
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA

nd_area <- dplyr::bind_rows(spec_preds)
area_res <- dplyr::bind_rows(res_list)

# Plot residuals
dfs <- data.frame(area = c("BS", "BT", "FB", "FM", "MU", "SI_EK", "RA", "JM", "HO", "SI_HA"),
                  df = c(rep("\u03BD = 6", 7), rep("\u03BD = 10", 3)))
  
area_res %>% 
  ggplot(aes(sample = res)) +
  stat_qq(size = 0.75, shape = 21, fill = NA) +
  facet_wrap(~area) +
  stat_qq_line() +
  labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
  geom_text(data = dfs, aes(label = df, hjust = -0.1, vjust = 1.25), inherit.aes = FALSE,
              x = -Inf, y = Inf, size = 2.6, color = "grey20") +
  theme(aspect.ratio = 1)

ggsave(paste0(home, "/figures/supp/qq_temp_area.pdf"), width = 17, height = 17, units = "cm", device = cairo_pdf)

# In order to have a lower temperature in before-nuclear times (without any data to inform that), we can use the nearby areas.. so FM informs BT prior to nuclear
nd_area_sub <- nd_area %>% 
  mutate(keep = "N",
         keep = ifelse(area == "FM" & year <= 1980, "Y", keep), # use FM instead of BT
         keep = ifelse(area == "SI_EK" & year <= 1972, "Y", keep)) %>% # use SI_EK instead of SI_HA
  filter(keep == "Y") %>% # Now change the labels to BT and SI_EK...
  mutate(area = ifelse(area == "FM", "BT", "SI_HA"))
mutate: new variable 'keep' (character) with 2 unique values and 0% NA
filter: removed 734,394 rows (90%), 81,252 rows remaining
mutate: converted 'area' from factor to character (0 new NA)
# Bind rows and plot only the temperature series we will use for growth modelling
nd_area <- bind_rows(nd_area, nd_area_sub) %>%
  select(-keep) %>% 
  mutate(growth_dat = ifelse(area == "SI_HA" & year %in% c(1966, 1967), "Y", growth_dat)) # SI_EK and SI_HA do not have the same starting years, so we can't use allo columns from SI_EK
select: dropped one variable (keep)
mutate: changed 2,196 values (<1%) of 'growth_dat' (0 new NA)
nd_area %>% 
  filter(area %in% c("FM", "BT", "SI_EK", "SI_HA")) %>% 
  filter(year <= 1980 & year >= 1966) %>% 
  group_by(area, year) %>% 
  summarise(mean_temp = mean(pred)) %>% 
  ungroup() %>% 
  pivot_wider(names_from = area, values_from = mean_temp)
filter: removed 532,362 rows (59%), 364,536 rows remaining
filter: removed 298,656 rows (82%), 65,880 rows remaining
group_by: 2 grouping variables (area, year)
summarise: now 60 rows and 3 columns, one group variable remaining (area)
ungroup: no grouping variables
pivot_wider: reorganized (area, mean_temp) into (BT, FM, SI_EK, SI_HA) [was 60x3, now 15x5]
# A tibble: 15 × 5
    year    BT    FM SI_EK SI_HA
   <dbl> <dbl> <dbl> <dbl> <dbl>
 1  1966  8.15  8.15  7.97  7.97
 2  1967  8.98  8.98  8.76  8.76
 3  1968  8.69  8.69  8.37  8.37
 4  1969  8.55  8.55  8.53  8.53
 5  1970  8.30  8.30  8.04  8.04
 6  1971  8.33  8.33  8.33  8.33
 7  1972  8.83  8.83  8.90  8.90
 8  1973  9.19  9.19  8.99 11.1 
 9  1974  9.05  9.05  8.66  8.46
10  1975  7.43  7.43  9.34 14.6 
11  1976  7.09  7.09  8.00 13.7 
12  1977  5.61  5.61  7.80 13.0 
13  1978  5.59  5.59  7.57 13.1 
14  1979  5.86  5.86  7.42 11.4 
15  1980  6.96  6.96  7.73 13.4 

Plot detailed exploration of predictions

# Loop trough all areas, plot temperature as a function of yday, color by data source, facet by year

nd_area %>% 
  ggplot(aes(yday, y = year, fill = pred, color = pred)) +
  geom_raster() +
  facet_wrap(~area, ncol = 4) +
  scale_fill_viridis(option = "magma") +
  scale_color_viridis(option = "magma") +
  labs(x = "Day-of-the-year", y = "Year", color = "Predicted SST (°C)", fill = "Predicted SST (°C)") + 
  theme(legend.position = c(0.78, 0.14))

ggsave(paste0(home, "/figures/supp/temp_pred_yday_area.pdf"), width = 17, height = 17, units = "cm")
ggsave(paste0(home, "/figures/for-talks/temp_pred_yday_area.pdf"), width = 14, height = 14, units = "cm")

for(i in unique(nd_area$area)) {
  
  plotdat <- nd_area %>% filter(area == i)
  
  print(
    ggplot(plotdat, aes(yday, pred, color = source)) + 
      scale_color_brewer(palette = "Dark2") + 
      facet_wrap(~year) + 
      geom_point(data = filter(d, area == i & year > min(plotdat$year)), size = 0.2,
                 aes(yday, temp, color = source)) + 
      geom_line(linewidth = 0.3) + 
      labs(title = paste("Site = ", i), color = "", linetype = "") + 
      guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
      theme_sleek(base_size = 8) +
      theme(legend.position = "bottom", 
            legend.direction = "horizontal",
            legend.spacing.y = unit(-0.3, "cm")) + 
      labs(x = "Day of the year", y = "Predicted SST (°C)")
  )
  
  ggsave(paste0(home, "/figures/supp/temp_pred_yday_area_", i, ".pdf" ), width = 17, height = 17, units = "cm")
  
}
filter: removed 813,234 rows (91%), 83,664 rows remaining
filter: removed 88,476 rows (96%), 3,449 rows remaining
filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 82,953 rows (90%), 8,972 rows remaining

filter: removed 806,013 rows (90%), 90,885 rows remaining
filter: removed 82,464 rows (90%), 9,461 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 79,969 rows (87%), 11,956 rows remaining

filter: removed 806,013 rows (90%), 90,885 rows remaining
filter: removed 83,775 rows (91%), 8,150 rows remaining

filter: removed 812,238 rows (91%), 84,660 rows remaining
filter: removed 81,803 rows (89%), 10,122 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 81,052 rows (88%), 10,873 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 87,117 rows (95%), 4,808 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 83,064 rows (90%), 8,861 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 76,748 rows (83%), 15,177 rows remaining

# Same but trimmed
for(i in unique(nd_area$area)) {
  
  plotdat <- nd_area %>% filter(area == i) %>% filter(growth_dat == "Y")
  
  print(
    ggplot(plotdat, aes(yday, pred, color = source)) + 
      scale_color_brewer(palette = "Dark2") + 
      facet_wrap(~year) + 
      geom_point(data = filter(d, area == i & year > min(plotdat$year)), size = 0.2,
                 aes(yday, temp, color = source)) + 
      geom_line(linewidth = 0.3) + 
      labs(title = paste("Site = ", i), color = "", linetype = "") + 
      guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
      theme_sleek(base_size = 8) +
      theme(#legend.position = c(0.8, 0.08), 
            legend.position = "bottom", 
            legend.direction = "horizontal",
            legend.spacing.y = unit(-0.3, "cm")) + 
      labs(x = "Day of the year", y = "Predicted SST (°C)")
  )
  
  ggsave(paste0(home, "/figures/supp/temp_pred_yday_area_trimmed_", i, ".pdf" ), width = 17, height = 17, units = "cm")
  
}
filter: removed 813,234 rows (91%), 83,664 rows remaining
filter: removed 65,520 rows (78%), 18,144 rows remaining
filter: removed 89,016 rows (97%), 2,909 rows remaining
filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 49,410 rows (54%), 41,724 rows remaining
filter: removed 82,953 rows (90%), 8,972 rows remaining

filter: removed 806,013 rows (90%), 90,885 rows remaining
filter: removed 38,325 rows (42%), 52,560 rows remaining
filter: removed 82,812 rows (90%), 9,113 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 49,410 rows (54%), 41,724 rows remaining
filter: removed 80,245 rows (87%), 11,680 rows remaining

filter: removed 806,013 rows (90%), 90,885 rows remaining
filter: removed 52,560 rows (58%), 38,325 rows remaining
filter: removed 84,279 rows (92%), 7,646 rows remaining

filter: removed 812,238 rows (91%), 84,660 rows remaining
filter: removed 19,380 rows (23%), 65,280 rows remaining
filter: removed 81,959 rows (89%), 9,966 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 69,174 rows (76%), 21,960 rows remaining
filter: removed 81,544 rows (89%), 10,381 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 66,978 rows (73%), 24,156 rows remaining
filter: removed 87,657 rows (95%), 4,268 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 38,430 rows (42%), 52,704 rows remaining
filter: removed 83,400 rows (91%), 8,525 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 37,332 rows (41%), 53,802 rows remaining
filter: removed 76,748 rows (83%), 15,177 rows remaining

Plot summarized data and predictions

# Summarise data
dsum <- d %>% 
  group_by(year, area, source) %>% 
  summarise(temp = mean(temp)) %>% 
  mutate(type = "data")
group_by: 3 grouping variables (year, area, source)
summarise: now 1,376 rows and 4 columns, 2 group variables remaining (year, area)
mutate (grouped): new variable 'type' (character) with one unique value and 0% NA
preds_area <- nd_area %>% 
  filter(growth_dat == "Y" & source == "logger") %>% 
  group_by(area, year) %>% 
  summarise(temp = mean(pred)) %>% 
  mutate(model = "area model")
filter: removed 760,105 rows (85%), 136,793 rows remaining
group_by: 2 grouping variables (area, year)
summarise: now 380 rows and 3 columns, one group variable remaining (area)
mutate (grouped): new variable 'model' (character) with one unique value and 0% NA

Make final plot using the area-specific model

# Overall t_opt from 02-fit-vbge.qmd! Update when I got final model!
topt_lwr <- 8.771866
topt <- 9.468607
topt_upr <- 10.304697

# Add latitude
area <- c("BS", "BT", "FB", "FM", "HO", "JM", "MU", "RA", "SI_EK", "SI_HA")
nareas <- length(area)
lat <- c(60, 60.4, 60.3, 60.5, 63.7, 58, 59, 65.9, 57.3, 57.4)
lon <- c(21.5, 18.1, 19.5, 18, 20.9, 16.8, 18.1, 22.3, 16.6, 16.7)
area_attr <- data.frame(cbind(area = area, lat = lat, lon = lon)) %>%
  mutate_at(c("lat","lon"), as.numeric) %>% 
  arrange(desc(lat))
mutate_at: converted 'lat' from character to double (0 new NA)
           converted 'lon' from character to double (0 new NA)
ggplot(preds_area, aes(year, temp, color = temp)) + 
  facet_wrap(~factor(area, levels = area_attr$area), ncol = 4) + 
  geom_hline(yintercept = topt_lwr, linewidth = 0.3, linetype = 2, color = "tomato3", alpha = 0.3) +
  geom_hline(yintercept = topt, linewidth = 0.3, linetype = 1, color = "tomato3", alpha = 0.3) +
  geom_hline(yintercept = topt_upr, linewidth = 0.3, linetype = 2, color = "tomato3", alpha = 0.3) +
  geom_line(linewidth = 0.6) +
  labs(x = "Year", y = "Model-predicted annual average temperature") + 
  scale_color_viridis(option = "magma", name = "Site") +
  guides(color = "none") 

ggsave(paste0(home, "/figures/annual_average_temperature.pdf"), width = 17, height = 12, units = "cm")
#ggsave(paste0(home, "/figures/for-talks/annual_average_temperature.pdf"), width = 15, height = 10, units = "cm")

# Check year range
preds_area %>% 
  group_by(area) %>% 
  summarise(min = min(year),
            max = max(year)) %>% 
  arrange(min)
group_by: one grouping variable (area)
summarise: now 10 rows and 3 columns, ungrouped
# A tibble: 10 × 3
   area    min   max
   <chr> <dbl> <dbl>
 1 JM     1953  2016
 2 BT     1963  2000
 3 FM     1963  2000
 4 SI_HA  1966  2014
 5 SI_EK  1968  2015
 6 FB     1969  2016
 7 MU     1981  2000
 8 HO     1982  2016
 9 BS     1985  2002
10 RA     1985  2006
# Check temperature range
preds_area %>% 
  group_by(area) %>% 
  summarise(min = min(temp),
            max = max(temp)) %>% 
  arrange(min)
group_by: one grouping variable (area)
summarise: now 10 rows and 3 columns, ungrouped
# A tibble: 10 × 3
   area    min   max
   <chr> <dbl> <dbl>
 1 RA     3.06  9.56
 2 JM     3.37 11.8 
 3 HO     3.99  8.07
 4 FB     5.04 10.6 
 5 MU     5.16  8.57
 6 SI_EK  5.49 10.4 
 7 BT     5.87 16.6 
 8 FM     5.87 15.0 
 9 BS     6.22 11.8 
10 SI_HA  7.76 16.9 
# Save prediction df
write_csv(preds_area, paste0(home, "/output/gam_predicted_temps.csv"))